<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of sppf</title>
  <meta name="keywords" content="sppf">
  <meta name="description" content="SPPF  Sigma-Point Particle Filter.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">core</a> &gt; sppf.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../menu.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\core&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>sppf
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>SPPF  Sigma-Point Particle Filter.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [estimate, ParticleFilterDS, pNoise, oNoise] = sppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> SPPF  Sigma-Point Particle Filter.

   This hybrid particle filter uses a sigma-point Kalman filter (SRUKF or SRCDKF) for proposal distribution generation
   and is an extension of the original &quot;Unscented Particle Filter&quot; of Van der Merwe, De Freitas &amp; Doucet.

   [estimate, ParticleFilterDS, pNoise, oNoise] = sppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)

   This filter assumes the following standard state-space model:

     x(k) = ffun[x(k-1),v(k-1),U1(k-1)]
     y(k) = hfun[x(k),n(k),U2(k)]

   where x is the system state, v the process noise, n the observation noise, u1 the exogenous input to the state
   transition function, u2 the exogenous input to the state observation function and y the noisy observation of the
   system.

   INPUT
         ParticleFilterDS     Particle filter data structure. (see field definitions below)
         pNoise               (NoiseDS) process noise data structure
         oNoise               (NoiseDS) observation noise data structure
         obs                  noisy observations starting at time k ( y(k),y(k+1),...,y(k+N-1) )
         U1                   exogenous input to state transition function starting at time k-1 ( u1(k-1),u1(k),...,u1(k+N-2) )
         U2                   exogenous input to state observation function starting at time k  ( u2(k),u2(k+1),...,u2(k+N-1) )
         InferenceDS          Inference data structure generated by GENINFDS function.

   OUTPUT
         estimate             State estimate generated from posterior distribution of state given all observation. Type of
                              estimate is specified by InferenceDS.estimateType
         ParticleFilterDS     Updated Particle filter data structure. Contains set of particles as well as their corresponding
                              weights.
         pNoise               process noise data structure     (possibly updated)
         oNoise               observation noise data structure (possibly updated)

   ParticleFilterDS fields:
         .N                   (scalar) number of particles
         .particles           (statedim-by-N matrix) particle mean buffer
         .particlesCov        (statedim-by-statedim-by-N matrix) particle covariance buffer (Cholesky factors)
         .pNoise              (NoiseDS) process noise data structure for SPKF based proposal generation
         .oNoise              (NoiseDS) observation noise data structure for SPKF based proposal generation
         .weights             (1-by-N r-vector) particle weights

   Required InferenceDS fields:
         .spkfType            (string) Type of SPKF to use (srukf or srcdkf).
         .estimateType        (string) Estimate type : 'mean', 'mode', etc.
         .resampleThreshold   (scalar) If the ratio of the 'effective particle set size' to the total number of particles
                                       drop below this threshold  i.e.  (N_efective/N) &lt; resampleThreshold
                                       the particles will be resampled.  (N_efective is always less than or equal to N)

   NOTE : All covariances are assumed to be of type 'sqrt', i.e. Cholesky factors.

   See also
   <a href="pf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = pf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">PF</a>, <a href="srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">SRUKF</a>, <a href="srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">SRCDKF</a>
   Copyright (c) Oregon Health &amp; Science University (2006)

   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
   academic use only (see included license file) and can be obtained from
   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
   software should contact rebel@csee.ogi.edu for commercial licensing information.

   See LICENSE (which should be part of the main toolkit distribution) for more
   detail.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>	CVECREP  Column vector replicate</li><li><a href="residualresample.html" class="code" title="function outIndex = residualResample(inIndex, weights);">residualresample</a>	RESIDUALRESAMPLE  Residual resampling for SIR. Performs the resampling stage of</li><li><a href="srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>	SRCDKF  Square Root Central Difference Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srukf</a>	SRUKF  Square Root Unscented Kalman Filter (Sigma-Point Kalman Filter variant)</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="../.././ReBEL-0.2.7/examples/state_estimation/demse2.html" class="code" title="">demse2</a>	DEMSE2  Demonstrate state estimation on a simple scalar nonlinear (time variant) problem</li><li><a href="../.././ReBEL-0.2.7/examples/state_estimation/demse4.html" class="code" title="">demse4</a>	DEMSE4  Bearing Only Tracking Example</li><li><a href="../.././ReBEL-0.2.7/examples/state_estimation/demse5.html" class="code" title="">demse5</a>	DEMSE4  Bearing and Frequency Tracking Example</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [estimate, ParticleFilterDS, pNoise, oNoise] = sppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)</a>
0002 
0003 <span class="comment">% SPPF  Sigma-Point Particle Filter.</span>
0004 <span class="comment">%</span>
0005 <span class="comment">%   This hybrid particle filter uses a sigma-point Kalman filter (SRUKF or SRCDKF) for proposal distribution generation</span>
0006 <span class="comment">%   and is an extension of the original &quot;Unscented Particle Filter&quot; of Van der Merwe, De Freitas &amp; Doucet.</span>
0007 <span class="comment">%</span>
0008 <span class="comment">%   [estimate, ParticleFilterDS, pNoise, oNoise] = sppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)</span>
0009 <span class="comment">%</span>
0010 <span class="comment">%   This filter assumes the following standard state-space model:</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%     x(k) = ffun[x(k-1),v(k-1),U1(k-1)]</span>
0013 <span class="comment">%     y(k) = hfun[x(k),n(k),U2(k)]</span>
0014 <span class="comment">%</span>
0015 <span class="comment">%   where x is the system state, v the process noise, n the observation noise, u1 the exogenous input to the state</span>
0016 <span class="comment">%   transition function, u2 the exogenous input to the state observation function and y the noisy observation of the</span>
0017 <span class="comment">%   system.</span>
0018 <span class="comment">%</span>
0019 <span class="comment">%   INPUT</span>
0020 <span class="comment">%         ParticleFilterDS     Particle filter data structure. (see field definitions below)</span>
0021 <span class="comment">%         pNoise               (NoiseDS) process noise data structure</span>
0022 <span class="comment">%         oNoise               (NoiseDS) observation noise data structure</span>
0023 <span class="comment">%         obs                  noisy observations starting at time k ( y(k),y(k+1),...,y(k+N-1) )</span>
0024 <span class="comment">%         U1                   exogenous input to state transition function starting at time k-1 ( u1(k-1),u1(k),...,u1(k+N-2) )</span>
0025 <span class="comment">%         U2                   exogenous input to state observation function starting at time k  ( u2(k),u2(k+1),...,u2(k+N-1) )</span>
0026 <span class="comment">%         InferenceDS          Inference data structure generated by GENINFDS function.</span>
0027 <span class="comment">%</span>
0028 <span class="comment">%   OUTPUT</span>
0029 <span class="comment">%         estimate             State estimate generated from posterior distribution of state given all observation. Type of</span>
0030 <span class="comment">%                              estimate is specified by InferenceDS.estimateType</span>
0031 <span class="comment">%         ParticleFilterDS     Updated Particle filter data structure. Contains set of particles as well as their corresponding</span>
0032 <span class="comment">%                              weights.</span>
0033 <span class="comment">%         pNoise               process noise data structure     (possibly updated)</span>
0034 <span class="comment">%         oNoise               observation noise data structure (possibly updated)</span>
0035 <span class="comment">%</span>
0036 <span class="comment">%   ParticleFilterDS fields:</span>
0037 <span class="comment">%         .N                   (scalar) number of particles</span>
0038 <span class="comment">%         .particles           (statedim-by-N matrix) particle mean buffer</span>
0039 <span class="comment">%         .particlesCov        (statedim-by-statedim-by-N matrix) particle covariance buffer (Cholesky factors)</span>
0040 <span class="comment">%         .pNoise              (NoiseDS) process noise data structure for SPKF based proposal generation</span>
0041 <span class="comment">%         .oNoise              (NoiseDS) observation noise data structure for SPKF based proposal generation</span>
0042 <span class="comment">%         .weights             (1-by-N r-vector) particle weights</span>
0043 <span class="comment">%</span>
0044 <span class="comment">%   Required InferenceDS fields:</span>
0045 <span class="comment">%         .spkfType            (string) Type of SPKF to use (srukf or srcdkf).</span>
0046 <span class="comment">%         .estimateType        (string) Estimate type : 'mean', 'mode', etc.</span>
0047 <span class="comment">%         .resampleThreshold   (scalar) If the ratio of the 'effective particle set size' to the total number of particles</span>
0048 <span class="comment">%                                       drop below this threshold  i.e.  (N_efective/N) &lt; resampleThreshold</span>
0049 <span class="comment">%                                       the particles will be resampled.  (N_efective is always less than or equal to N)</span>
0050 <span class="comment">%</span>
0051 <span class="comment">%   NOTE : All covariances are assumed to be of type 'sqrt', i.e. Cholesky factors.</span>
0052 <span class="comment">%</span>
0053 <span class="comment">%   See also</span>
0054 <span class="comment">%   PF, SRUKF, SRCDKF</span>
0055 <span class="comment">%   Copyright (c) Oregon Health &amp; Science University (2006)</span>
0056 <span class="comment">%</span>
0057 <span class="comment">%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for</span>
0058 <span class="comment">%   academic use only (see included license file) and can be obtained from</span>
0059 <span class="comment">%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the</span>
0060 <span class="comment">%   software should contact rebel@csee.ogi.edu for commercial licensing information.</span>
0061 <span class="comment">%</span>
0062 <span class="comment">%   See LICENSE (which should be part of the main toolkit distribution) for more</span>
0063 <span class="comment">%   detail.</span>
0064 
0065 <span class="comment">%=============================================================================================</span>
0066 <span class="keyword">if</span> (nargin ~= 7) error(<span class="string">' [ sppf ] Not enough input arguments.'</span>); <span class="keyword">end</span>
0067 
0068 <span class="comment">%--------------------------------------------------------------------------------------------------</span>
0069 
0070 Xdim  = InferenceDS.statedim;                            <span class="comment">% extract state dimension</span>
0071 Odim  = InferenceDS.obsdim;                              <span class="comment">% extract observation dimension</span>
0072 U1dim = InferenceDS.U1dim;                               <span class="comment">% extract exogenous input 1 dimension</span>
0073 U2dim = InferenceDS.U2dim;                               <span class="comment">% extract exogenous input 2 dimension</span>
0074 Vdim  = InferenceDS.Vdim;                                <span class="comment">% extract process noise dimension</span>
0075 Ndim  = InferenceDS.Ndim;                                <span class="comment">% extract observation noise dimension</span>
0076 
0077 N  = ParticleFilterDS.N;                      <span class="comment">% number of particles</span>
0078 x  = ParticleFilterDS.particles;              <span class="comment">% copy particle mean buffer</span>
0079 Sx = ParticleFilterDS.particlesCov;           <span class="comment">% copy particle covariance buffer</span>
0080 pNoiseSPKF = ParticleFilterDS.pNoise;         <span class="comment">% copy SPKF process noise data structure</span>
0081 oNoiseSPKF = ParticleFilterDS.oNoise;         <span class="comment">% copy SPKF observation noise data structure</span>
0082 weights    = ParticleFilterDS.weights;        <span class="comment">% particle weights</span>
0083 St         = round(N*InferenceDS.resampleThreshold);   <span class="comment">% resample threshold</span>
0084 
0085 onoise = zeros(InferenceDS.Ndim,N);
0086 
0087 normWeights = <a href="cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(1/N,N);
0088 
0089 NOV = size(obs,2);                                       <span class="comment">% number of input vectors</span>
0090 
0091 <span class="keyword">if</span> (U1dim==0), UU1=zeros(0,N); <span class="keyword">end</span>
0092 <span class="keyword">if</span> (U2dim==0), UU2=zeros(0,N); <span class="keyword">end</span>
0093 
0094 estimate   = zeros(Xdim,NOV);
0095 
0096 SxPred  = zeros(Xdim,Xdim,N);                     <span class="comment">% setup buffers</span>
0097 xNew    = zeros(Xdim,N);
0098 xPred   = zeros(Xdim,N);
0099 
0100 ones_numP = ones(N,1);
0101 ones_Xdim = ones(1,Xdim);
0102 
0103 proposal = zeros(1,N);
0104 
0105 normfact = (2*pi)^(Xdim/2);
0106 
0107 
0108 <span class="keyword">for</span> j=1:NOV,
0109 <span class="comment">%---------------------------------------------------------------------------</span>
0110 
0111     OBStemp = obs(:,j);                <span class="comment">% inline cvecrep</span>
0112     OBS = OBStemp(:,ones_numP);
0113 
0114     <span class="keyword">if</span> U1dim
0115       Utemp = U1(:,j);
0116       UU1 = Utemp(:,ones_numP);        <span class="comment">% inline cvecrep</span>
0117     <span class="keyword">end</span>
0118     <span class="keyword">if</span> U2dim
0119       Utemp = U2(:,j);
0120       UU2 = Utemp(:,ones_numP);        <span class="comment">% inline cvecrep</span>
0121     <span class="keyword">end</span>
0122 
0123 
0124     <span class="comment">%-----------------------------------------------------------------------</span>
0125     <span class="comment">% TIME UPDATE</span>
0126 
0127     randBuf = randn(Xdim,N);
0128 
0129     <span class="keyword">switch</span> InferenceDS.spkfType
0130     <span class="keyword">case</span> <span class="string">'srukf'</span>
0131         <span class="keyword">for</span> k=1:N,
0132             [xNew(:,k), SxPred(:,:,k), pNoiseSPKF, oNoiseSPKF, intvarDS] = <a href="srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srukf</a>(x(:,k), Sx(:,:,k), pNoiseSPKF, oNoiseSPKF, <span class="keyword">...</span>
0133                                                                         obs(:,j), UU1(:,k), UU2(:,k), InferenceDS);
0134 
0135             xPred(:,k) = xNew(:,k) + SxPred(:,:,k)*randBuf(:,k);
0136         <span class="keyword">end</span>
0137     <span class="keyword">case</span> <span class="string">'srcdkf'</span>
0138         <span class="keyword">for</span> k=1:N,
0139             [xNew(:,k), SxPred(:,:,k), pNoiseSPKF, oNoiseSPKF, intvarDS] = <a href="srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>(x(:,k), Sx(:,:,k), pNoiseSPKF, oNoiseSPKF, <span class="keyword">...</span>
0140                                                                         obs(:,j), UU1(:,k), UU2(:,k), InferenceDS);
0141             xPred(:,k) = xNew(:,k) + SxPred(:,:,k)*randBuf(:,k);
0142         <span class="keyword">end</span>
0143     <span class="keyword">otherwise</span>
0144         error(<span class="string">' [ sppf ] Unknown SPKF type.'</span>);
0145     <span class="keyword">end</span>
0146 
0147     <span class="comment">%-----------------------------------------------------------------------</span>
0148     <span class="comment">% EVALUATE IMPORTANCE WEIGHTS</span>
0149 
0150 
0151     <span class="comment">% calculate transition prior for each particle (in log domain)</span>
0152     prior = InferenceDS.prior( InferenceDS, xPred, x, UU1, pNoise) + 1e-99;
0153 
0154     <span class="comment">% calculate observation likelihood for each particle (in log domain)</span>
0155     likelihood = InferenceDS.likelihood( InferenceDS, OBS, xPred, UU2, oNoise) + 1e-99;
0156 
0157     difX = xPred - xNew;
0158 
0159 
0160     <span class="keyword">for</span> k=1:N,
0161         cholFact = SxPred(:,:,k);
0162         foo = cholFact \ difX(:,k);
0163         proposal(k) = exp(-0.5*foo'*foo) / abs(normfact*prod(diag(cholFact))) + 1e-99;
0164         weights(k) = weights(k) * likelihood(k) * prior(k) / proposal(k);
0165     <span class="keyword">end</span>
0166 
0167 
0168 <span class="keyword">if</span> (0)
0169 figure(20);
0170 subplot(411)
0171 stem(prior);
0172 ylabel(<span class="string">'prior'</span>);
0173 subplot(412);
0174 stem(likelihood);
0175 ylabel(<span class="string">'likelihood'</span>);
0176 subplot(413);
0177 stem(proposal);
0178 ylabel(<span class="string">'proposal'</span>);
0179 subplot(414);
0180 stem(weights);
0181 ylabel(<span class="string">'weights'</span>);
0182 drawnow
0183 <span class="keyword">end</span>
0184 
0185     weights = weights / sum(weights);
0186 
0187     <span class="comment">%-----------------------------------------------------------------------</span>
0188     <span class="comment">% CALCULATE ESTIMATE</span>
0189 
0190     <span class="keyword">switch</span> InferenceDS.estimateType
0191 
0192     <span class="keyword">case</span> <span class="string">'mean'</span>
0193         muFoo = sum(weights(ones_Xdim,:).*xPred,2);
0194         estimate(:,j) = muFoo;          <span class="comment">% expected mean</span>
0195 
0196     <span class="keyword">otherwise</span>
0197         error(<span class="string">' [ sppf ] Unknown estimate type.'</span>);
0198 
0199     <span class="keyword">end</span>
0200 
0201 
0202     <span class="comment">%-----------------------------------------------------------------------</span>
0203     <span class="comment">% RESAMPLE</span>
0204 
0205     S = 1/sum(weights.^2);     <span class="comment">% calculate effective particle set size</span>
0206 
0207     <span class="keyword">if</span> (S &lt; St)                  <span class="comment">% resample if S is below threshold</span>
0208         outIndex  = <a href="residualresample.html" class="code" title="function outIndex = residualResample(inIndex, weights);">residualresample</a>(1:N,weights);
0209         x = xPred(:,outIndex);
0210         <span class="keyword">for</span> k=1:N,
0211             Sx(:,:,k) = SxPred(:,:,outIndex(k));
0212         <span class="keyword">end</span>
0213         weights = normWeights;
0214     <span class="keyword">else</span>
0215         x  = xPred;
0216         Sx = SxPred;
0217     <span class="keyword">end</span>
0218 
0219 
0220 
0221     <span class="keyword">if</span> pNoise.adaptMethod <span class="keyword">switch</span> InferenceDS.inftype
0222     <span class="comment">%---------------------- UPDATE PROCESS NOISE SOURCE IF NEEDED --------------------------------------------</span>
0223     <span class="keyword">case</span> <span class="string">'parameter'</span>  <span class="comment">%--- parameter estimation</span>
0224         <span class="keyword">switch</span> pNoise.adaptMethod
0225         <span class="keyword">case</span> <span class="string">'anneal'</span>
0226             pNoise.cov = diag(max(pNoise.adaptParams(1) * diag(pNoise.cov) , pNoise.adaptParams(2)));
0227         <span class="keyword">otherwise</span>
0228             error(<span class="string">' [ sppf ]unknown process noise adaptation method!'</span>);
0229         <span class="keyword">end</span>
0230 
0231     <span class="keyword">case</span> <span class="string">'joint'</span>  <span class="comment">%--- joint estimation</span>
0232         idx = pNoise.idxArr(<span class="keyword">end</span>,:); <span class="comment">% get indexs of parameter block of combo noise source</span>
0233         idxRange=idx(1):idx(2);
0234         <span class="keyword">switch</span> pNoise.adaptMethod
0235         <span class="keyword">case</span> <span class="string">'anneal'</span>
0236             pNoise.cov(idxRange,idxRange) = diag(max(pNoise.adaptParams(1) * diag(pNoise.cov(idxRange,idxRange)), pNoise.adaptParams(2)));
0237         <span class="keyword">otherwise</span>
0238             error(<span class="string">' [ sppf ]unknown process noise adaptation method!'</span>);
0239         <span class="keyword">end</span>
0240         pNoise.noiseSources{pNoise.N}.cov = pNoise.cov(idxRange,idxRange);
0241 
0242     <span class="comment">%--------------------------------------------------------------------------------------------------</span>
0243     <span class="keyword">end</span>; <span class="keyword">end</span>
0244 
0245 <span class="comment">%--------------------------------------------------------------------------</span>
0246 <span class="keyword">end</span>     <span class="comment">%.. loop over input vectors</span>
0247 
0248 
0249 ParticleFilterDS.particles    = x;
0250 ParticleFilterDS.particlesCov = Sx;
0251 ParticleFilterDS.weights      = weights;
0252 ParticleFilterDS.pNoise       = pNoiseSPKF;
0253 ParticleFilterDS.oNoise       = oNoiseSPKF;</pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>